
%\documentclass{cumcmthesis}
\documentclass[withoutpreface,bwprint]{cumcmthesis} %去掉封面与编号页
\usepackage{siunitx} %用于摄氏度
\usepackage{graphicx} %图片需要使用的宏包
\usepackage{subfigure}	%用于排版多张图片
\usepackage{float}	%用于排版图片位置
\bibliographystyle{plain}	%引用样式,参考文献

\usepackage{url}
\title{基于一维热传导方程的炉温曲线优化问题}

\begin{document}

\maketitle
\begin{abstract}
	在集成电路板等电子产品生产中，需要将安装有各种电子元件的印刷电路板放置在回焊炉中，通过加热，将电子元件自动焊接到电路板上。在这个生产过程中，让回焊炉的各部分保持工艺要求的温度，对产品质量至关重要。目前，这方面的许多工作是通过实验测试来进行控制和调整的。本文旨在通过机理模型来进行分析研究。
	\par
	对于问题一，使用一维热传导方程求解炉内各个小温区间温度，得出炉内温度分布。然后，对于印刷电路板温度变化过程，再利用一维热传导方程，以牛顿冷却定律修正后的焊接区域表面温度为边界条件，拟合函数曲线，求解出热传导方程、冷却定律中的最优参数。最后得到焊接区域中心温度变化规律，并做出炉温曲线，将数据存入result.csv中，得到的小温区3、6、7中点及小温区8结束处焊接区域中心的温度分别为128.98℃、166.65 ℃、188.98℃、223.28℃。
	\par
	对于问题二，为求出传送带的最大速度，我们先根据炉内温度分布模型将问题二中的炉内温度分布求解出来。我们先根据炉内温度分布模型将问题二中的炉内温度分布求解出来。在制程界限范围内，以最大传送带速度为目标函数，建立单目标优化模型，使用遍历方法，求解出最大速度为80.8cm/min。
	\par
	对于问题三，为了使阴影面积最小，我们先建立阴影面积求解公式，利用正弦曲线拟合温度曲线，简化阴影面积求解公式。然后以最小阴影面积为目标函数、制程界限为约束条件的优化模型。最后通过模拟退火算法求解出最优参数(185,185,235,265,95)，最小阴影面积为525.44，并得到最优炉温曲线。
	\par
	对于问题四，在第三问的基础上，为了使图形对称程度较高，我们建立两个指标，分别判断阴影面积的大小、曲线的对称程度，将阴影面积大小作为约束条件，曲线尽可能对称作为目标函数，建立单目标优化模型，使用模拟退火算法求解出最优参数(185,193,230,265,94)。
	
\keywords{ \textbf{热传导方程}\quad  \textbf{牛顿冷却定律} \quad \textbf{优化模型} \quad \textbf{模拟退火} }

\end{abstract}

%目录
%\tableofcontents

%新的一页
%\newpage

\section{问题背景}

\subsection{问题背景}

在生产集成电路板等电子产品的过程中，需要将各种电子元件的印刷电路板在回焊炉内进行焊接。通过加热的方式将电子元件自动焊接到电路板上。在此过程中对回焊炉内的温度、电路板两侧传送带的速度和电子元件的焊接时间等因素有着严格要求，对产品质量至关重要。实际生产时可以通过调节各温区的设定温度和传送带的过炉速度来控制产品质量。目前，这方面的许多工作是通过实验测试来进行控制和调整的。通常情况下，回焊炉内部设置若干个小温区，它们从功能上可分成４个大温区：预热区、恒温区、回流区、冷却区．电路板两侧搭在传送带上匀速进入炉内进行加热焊接．某回焊炉内有11个小温区及炉前区域和炉后区域。设定温度和传送带速度控制在一定的范围内并且炉温曲线需要满足制程界限。



\section{问题重述}

\subsection{问题一}

假设传送带过炉速度为78 cm/min，各温区温度的设定值分别为173ºC（小温区1~5）、198ºC（小温区6）、230ºC（小温区7）和257ºC（小温区8~9），请对焊接区域的温度变化规律建立数学模型，列出小温区3、6、7中点及小温区8结束处焊接区域中心的温度，画出相应的炉温曲线，并将每隔0.5s焊接区域中心的温度存放在提供的result.csv中。

\subsection{问题二}

假设各温区温度的设定值分别为182ºC（小温区1-5）、203ºC（小温区6）、237ºC（小温区7）、254ºC（小温区8-9），请确定允许的最大传送带过炉速度。

\subsection{问题三}

为使炉温曲线理想，应使超过217ºC到峰值温度所覆盖的面积最小。请确定最优炉温曲线，以及各温区的设定温度和传送带的过炉速度，并给出相应的面积。

\subsection{问题四}

在焊接过程中，希望以峰值温度为中心线的两侧超过217ºC的炉温曲线应尽量对称。请结合问题3，进一步给出最优炉温曲线，以及各温区设定的温度及传送带过炉速度，并给出相应的指标值。


\section{模型的假设}

本文提出以下合理假设：

\begin{itemize}
% \item 假设该线路上运行的是同一钟类型的公交车；
\item 每个小温区内及边界都按设定的温度计算；
\item 印刷电路板在传送带上运动不会影响回焊炉内温度分布；
\item 每两个小温区之间的间隙的温度分布只收相邻两小温度区温度影响；

\end{itemize}

% 讲解三线表以及网站使用

% 讲解符号查询以及对应含义

\section{符号说明}
\begin{center}
\begin{tabular}{cc}
 \toprule[1.5pt]
 \makebox[0.3\textwidth][c]{符号}	&  \makebox[0.4\textwidth][c]{意义} \\
 \midrule[1pt]
%  $ W $	    	& 某一小时内该路段运行总收益-总成本   \\ 
%  $ W_0 $	    & 区分高峰和低峰的一个临界值  \\ 
%  $ P $	    	& 线路在一小时内所有站的总上车人数 \\ 
%  $ x $	    	& 线路在一小时内的车辆数 \\  
%  $ T_t $	    & 长期趋势项 \\ 
%  $ M_t $	    & 简单移动平均项 \\ 
    $T_1$ & \text{小温区 $1$-$5$ 的温度}\\
	$T_2$ & \text{小温区 $6$ 的温度}\\
	$T_3$ & \text{小温区 $7$ 的温度}\\
	$T_4$ & \text{小温区 $8$-$9$ 的温度}\\
	$T_5$ & \text{小温区 $10$-$11$ 的温度}\\
	$u_{i,k}$ &\text{印刷电路板焊接区域中心在第$i$个位置第$k$时的温度}\\
	$v$   & \text{传送带过炉速度}\\
	$u_{stove}$   & \text{炉内温度}\\
	$k$   & \text{牛顿冷却定律系数}\\
	$Q_{\eta }$& \text{阴影面积}   \\
\bottomrule[1.5pt]
\end{tabular}
\end{center}

\section{问题分析}

\subsection{问题一}
要对焊接区域的温度变化规律建立数学模型，并对问题求解。首先要将炉内各个小温区间间隙温度分布通过简化的一维热传导方程求解出来，得到炉内温度分布。然后，对于印刷电路板温度变化过程，再利用一维热传导方程，以牛顿冷却定律修正后的焊接区域表面温度为边界条件，通过拟合函数曲线求解出最优的参数。最后得到焊接区域中心温度变化规律，并做出炉温曲线，将数据存入result.csv

\subsection{问题二}
在各个小温区温度都给定的情况下，求解出最大的传送带过炉速度v。我们先根据炉内温度分布模型将问题二中的炉内温度分布求解出来。并根据问题一中的约束条件，我们建立出以传送带速度最大为目标函数的优化模型。

\subsection{问题三}
需要在满足制程界限的条件下使得理想的炉温曲线应使超过217ºC到峰值温度所覆盖的面积最小。我们先利用正弦曲线来近似模拟原曲线，进一步简化得到阴影面积的求解公式。然后以最小阴影面积为目标函数、制程界限为约束条件的优化模型。最后通过模拟退火算法求解出最小阴影面积，并得到最优炉温曲线。

\subsection{问题四}
在第三问的基础上，又附加一个为使阴影面积尽量对称。通过定义两个指标分别判断阴影面积的大小、曲线的对称程度，将阴影面积的大小作为约束条件，将问题转化成单目标优化问题，再通过模拟退火算法求解。


\section{问题一的模型建立及求解}

\subsection{模型建立}

\subsubsection{炉内温度分布模型}
根据假设，每个小温区内及边界都按设定的温度计算，我们只需要建立炉内各个小温区间间隙温度分布就能得到炉内温度分布模型。
根据热传导的傅里叶定律及能量守恒定律易推 出三维热传导方程

% 讲解公式以及引用
\begin{equation}
	\begin{gathered}
		\frac{\partial u}{\partial t}=a^{2}\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}+\frac{\partial^{2} u}{\partial z^{2}}\right) 
	\end{gathered}
\end{equation}

\begin{equation}
	a^{2}=\frac{k}{C \rho}
\end{equation}
其中，$k$，$C$，$\rho $分别为导热介质的热传导率、比热、密度。
对于炉内各个小温区间间隙温度分布，我们将其简化成沿传送带方向的一维热传导
问题。
由上式可得一维热传导方程为
\begin{equation}
	\frac{\partial u}{\partial t}=a^{2} \frac{\partial^{2} u}{\partial x^{2}}
\end{equation}
根据题目已知，焊接工作是在回焊炉内空气温度达到稳定后进行，此时炉内温度不随时间变化，炉内温度对时间的偏导为$0$，及
\begin{equation}
	\frac{\partial u}{\partial t}=0
\end{equation}

所以，
\begin{equation}
	\frac{\partial^2u}{\partial x^2}=0	
\end{equation}
求解得到$u = k x + b$。

即各小温区间间隙温度分布成线性分布。根据题目给出的回焊炉的尺寸，在已知各小温区温度的前提下，由线性模型得出回焊炉内温度$u_{stove}$分布曲线

\begin{figure}[H]
	\centering
	\includegraphics[width=14cm]{炉内温度分布.png} %[图片大小]{图片路径}
	\caption{炉内温度曲线} %图片标题
\end{figure}

\subsubsection{印刷电路板加热、冷却过程模型}
根据题目，已知印刷电路板通过传送带进入回焊炉，经过预热区、恒温区、回流区、冷却区。在这个过程中，印刷电路板经过先升温后降温这样一个过程。我们将印刷电路板温度变化过程看成以焊接区域表面一点到其焊接区域中心这样一个一维的热传到问题。
由1中推导出的一维热传导方程
\begin{equation}
	\frac{\partial u}{\partial t}=a^2\frac{\partial^2u}{\partial x^2}
\end{equation}


边界条件：
我们引入牛顿冷却定律
\begin{equation}
	\frac{du}{dt}=k(u_{stove}-u)	
\end{equation}

修正印刷电路板焊接区域边界条件
\par
初始条件：

\begin{equation}
	u=25
\end{equation}

\subsection{问题求解}

我们将$T_1$= \SI{175}{\degreeCelsius} 、$T_2$= \SI{195}{\degreeCelsius} 、$T_3$= \SI{235}{\degreeCelsius} 、$T_4$= \SI{255}{\degreeCelsius} 带入炉内温度曲线，我们利用原始数据求解出炉内温度曲线，对于推导出来的:                                           
\begin{equation}
	\frac{\partial u}{\partial t}=a^2\frac{\partial^2u}{\partial x^2}
\end{equation}
由于建立的热传导方程属于抛物型方程，边值条件复杂难以求得解析解。采
用有限差分法，将连续的定解区域用有效个离散点构成的网络来代替，把定解区
域上的连续函数用网格上定义的离散变量的函数来近似，把原方程和定解条件中
的微商用差商来近似。通过迭代求得参数下的电路板温度，最后通过借鉴最小二乘法
的思想确定最优参数。
\begin{equation}
	\begin{aligned}
		\begin{split}
		\left \{
		\begin{array}{ll}
		   \frac{\partial u(x,t)}{\partial t} |_{i,k} = \frac{u_{i,k+1}-u_{i,k}}{\tau } \\
			\\
		   \left.\frac{\partial^{2} u(x, t)}{\partial x^{2}}\right|_{i, k}=\frac{u_{i-1, k}-2 u_{i, k}+u_{i+1, k}}{h^{2}}\\
		\end{array}
		\right.
		\end{split}
	\end{aligned}
\end{equation}
得到热传导模型的差分方程为：
\begin{equation}
	u_{i, k+1}-u_{i, k}=\frac{\tau}{h^{2}} \cdot a^2 \cdot\left(u_{i-1, k}-2   u_{i, k}+u_{i+1, k}\right) 
\end{equation}
边界条件：
\begin{equation}
	\left \{
	\begin{array}{l}
  	u_{i, 0}=25 \\
  	\\
	u_{0,k+1}=-k \cdot \tau  \cdot(u_{0,k}-tem(k\tau v))+u_{0,k}\\
	\\
	u_{0.15,k+1}=-k \cdot \tau  \cdot(u_{0.15,k}-tem(k\tau v))+u_{0.15,k}\\
  	% u_{0, k}=\operatorname{tem}(k, v)
	\end{array}
	\right.
\end{equation}
其中$k\tau v$表示当前印刷电路板在回焊炉所处的位置，$tem(k\tau v)$表示印刷电路板所处位置环境的温度

我们将印刷电路板在回焊炉中温度变化过程分为五段，分别对应参数$a^2_1$，$a^2_2$，$a^2_3$，$a^2_4$，$a^2_5$，和由牛顿冷却定律修正后的边界条件中参数$k$。我们首先确定这六个参数的大致范围，然后在各自的范围内枚举这六个参数，得到最优值
\begin{table}[H]
	\centering
	\caption{参数值表}
	\begin{tabular}{ccc}
		\toprule[1.5pt]
		\makebox[0.3\textwidth][c]{参数}	& \makebox[0.4\textwidth][c]{数值范围} &  \makebox[0.4\textwidth][c]{数值} \\ 
		\midrule[1pt]		
		$k$   & 3000-6000 & \text{5000}\\
    	$a^2_1$ & 0.00003-0.00004 & \text{0.000038}\\
		$a^2_2$ & 0.00004-0.00005 & \text{0.000045}\\
		$a^2_3$ & 0.00006-0.00007 & \text{0.000063}\\
		$a^2_4$ & 0.00004-0.00005 & \text{0.000047}\\
		$a^2_5$ & 0.00002-0.00003 & \text{0.000022}\\
		\bottomrule[1.5pt]		
	\end{tabular}
\end{table}

通过遍历 $a$，得到$t$对应的温度$u$,然后与附件中的对应的数据求平方差和，以此作为决策依据，其中使得决策依据最小的我们认为是最优的参数$a$，并以这个参数绘出如下曲线，由此可以看出我们找到的$a$，拟合效果非常好。

\begin{figure}[H]
\centering
\includegraphics[width=17cm]{1.png} %[图片大小]{图片路径}
\caption{问题一拟合曲线} %图片标题
\end{figure}
将问题一中的温度$T_1$= \SI{173}{\degreeCelsius} 、$T_2$= \SI{198}{\degreeCelsius} 、$T_3$= \SI{230}{\degreeCelsius} 、$T_4$= \SI{257}{\degreeCelsius} 代入模型中，求解出炉温曲线如下
\begin{figure}[H]
	\centering
	\includegraphics[width=18cm]{第一问曲线.png} %[图片大小]{图片路径}
	\caption{问题一曲线} %图片标题
\end{figure}


并且，我们得到的曲线满足制程界限:

\begin{table}[H]
	\centering
	\caption{参数值表}
	\begin{tabular}{ccc}
		\toprule[1.5pt]
		\makebox[0.3\textwidth][c]{参数}	&  \makebox[0.4\textwidth][c]{数值} \\ 
		\midrule[1pt]		
		名称                & 值      & 单位    \\ 
		\hline
		最大温度上升斜率         & 2.002  & \SI{}{\degreeCelsius}/s \\
		\hline
		最小温度上升斜率           & -1.770 & \SI{}{\degreeCelsius}/s \\
		\hline
		温度上升过程中在150\SI{}{\degreeCelsius}-190\SI{}{\degreeCelsius}的时间 & 75.5   &    s    \\
		\hline
		温度大于\SI{217}{\degreeCelsius}的时间	 &   67   & s  \\
		\hline
		峰值温度	   &   240.67   &      \SI{}{\degreeCelsius}\\
		\bottomrule[1.5pt]		
	\end{tabular}
\end{table}



则小温区3、6、7中点及小温区8结束处焊接区域中心的温度分别为

\begin{center}
	\begin{tabular}{|c|c|}
		\hline \multicolumn{1}{|c|} { 位置 } & 温度 \\
		\hline 小温区 3 中点焊接区域中心 & $128.98^{\circ} \mathrm{C}$ \\
		\hline 小温区 6 中点焊接区域中心 & $166.65^{\circ} \mathrm{C}$ \\
		\hline 小温区 7 中点焊接区域中心 & $188.98^{\circ} \mathrm{C}$ \\
		\hline 小温区 8 结束处焊接区域中心 & $223.28^{\circ} \mathrm{C}$ \\
		\hline
	\end{tabular}
\end{center}

% \begin{tabular}{|c|c|c|}
% 	\hline 2&9\\
% 	\hline 7&5\\
% 	\hline 6&1\\
% 	\hline
% \end{tabular}
\section{问题二的模型建立及求解}
\subsection{模型建立}
该问题是在满足约束条件的前提下，尽可能的提高传送带的速度，我们建立一个单目标优化模型。
\subsubsection{优化模型}
$\bullet$ 决策变量
\begin{equation}
	\text{传送带速度}\quad v
\end{equation}

$\bullet$目标函数
\begin{equation}
	Max\quad v
\end{equation}

$\bullet$约束条件
\par
根据制程界限得出如下约束条件
\begin{equation}
	\begin{aligned}
		\begin{split}
		\left \{
		\begin{array}{ll}
			\frac{du}{dt}\in[-3,3]\\
			60 \leq\left. t\right|_{150 \leq u \leq 190} \leq 120 \\
			40 \leq\left. t\right|_{u \geq 270} \leq 90\\
			% when \quad u\in[150,190],t\in[60,120]\\
			% when \quad u\geq217,t\in[40,90]\\
			u_{max}\in[240,250]\\
		\end{array}
		\right.
		\end{split}
	\end{aligned}
\end{equation}

\subsection{问题求解}
速度v的制程界限为65 - 100，我们通过遍历，进行求解，最终求得最大速度为80.8cm/s
\begin{figure}[H]
	\centering
	\includegraphics[width=16cm]{第二问温度曲线图.png} %[图片大小]{图片路径}
	\caption{问题二曲线} %图片标题
\end{figure}


并且满足制程界限
\begin{table}[H]
	\centering
	\caption{制程界限}
	\begin{tabular}{ccc}
		\toprule[1.5pt]
		\makebox[0.3\textwidth][c]{参数}	&  \makebox[0.4\textwidth][c]{数值} \\ 
		\midrule[1pt]		
		名称                & 值      & 单位    \\ 
		\hline
		最大温度上升斜率         & 2.002  & \SI{}{\degreeCelsius}/s \\
		\hline
		最小温度上升斜率           & -1.770 & \SI{}{\degreeCelsius}/s \\
		\hline
		温度上升过程中在150\SI{}{\degreeCelsius}-190\SI{}{\degreeCelsius}的时间 & 75.5   &    s    \\
		\hline
		温度大于\SI{217}{\degreeCelsius}的时间	 &   67   & s  \\
		\hline
		峰值温度	   &   240.67   &      \SI{}{\degreeCelsius}\\
		\bottomrule[1.5pt]		
	\end{tabular}
\end{table}
最大速度80.8 cm/min。


\section{问题三的模型建立及求解}
\subsection{模型建立}
\subsubsection{阴影面积}
如图所示的阴影面积$Q_\eta$:
\begin{figure}[H]
	\centering
	\includegraphics[width=12cm]{第三问阴影面积.png} %[图片大小]{图片路径}
	\caption{阴影面积} %图片标题
\end{figure}
\begin{equation}
	Q_\eta=\int_{t_1}^{t_2}{(u(t)-T_M})dt
\end{equation}
\subsubsection{优化模型}
$\bullet$ 决策变量
\begin{equation}
	T=(T_1,T_2,T_3,T_4,v)
\end{equation}
 

$\bullet$ 目标函数
\begin{equation}
min\quad{Q_\eta(T)}
\end{equation}

$\bullet$ 约束条件
\begin{equation}
	\begin{aligned}
		\begin{split}
		\left \{
		\begin{array}{ll}
		    T_1\in[165,185]\\
			T_2\in[185,205]\\
			T_3\in[225,245]\\
			T_4\in[245,265]\\
			v\in[65,100]\\
			\frac{dT}{dt}\in[-3,3]\\
			60 \leq\left. t\right|_{150 \leq u \leq 190} \leq 120 \\
			40 \leq\left. t\right|_{u \geq 270} \leq 90\\
			u_{max}\in[240,250]\\
		\end{array}
		\right.
		\end{split}
	\end{aligned}
\end{equation}



\subsubsection{问题求解}
通过查阅文献，我们建立一个从$t_1-t_2$ ，高为 $\Delta T$ 的正半周期正弦曲线来近似模拟原曲线，将面积公式简化为：
\begin{equation}
	\begin{aligned}
		&Q_{\eta}=\int_{t_{1}}^{t_{2}} \Delta T \sin t d t \\
		&\quad=\frac{\Delta t \cdot \Delta T}{\pi} \int_{0}^{\pi} \sin t d t\\
		&\quad=\frac{2}{\pi} \Delta t \cdot \Delta T
	\end{aligned}
\end{equation}
\begin{figure}[H]
	\centering
	\includegraphics[width=12cm]{4.png} %[图片大小]{图片路径}
	\caption{示例模型} %图片标题
\end{figure}
模拟退火算法由N. Metropolis等人于1953年提出。1983年,S. Kirkpatrick等成功地将退火思想引入到组合优化领域。它是基于Monte-Carlo 迭代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合一定的概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。\par
其中最优解的更新机制为：如果新解比当前解更优，则接受新解，否则基于Metropolis准则判断是否接受新解。接受概率为：
\begin{equation}
	P= \begin{cases}1, & E_{t+1}<E_{t} \\ e^{\frac{-\left(E_{t+1}-E_{t}\right)}{k T}}, & E_{t+1} \geq E_{t}\end{cases}
\end{equation}
\par
如上公式：假设当前时刻搜索的解为$x_t$，对应的系统能量(目标函数)为$E_t$，对搜索点施加随机扰动，产生新解$x_{t+1}$，相应地，系统能量$E_{t+1}$，那么系统对搜索点从$x_t$到$x_{t+1}$转变的接受概率就为公式。
\\
使用模拟退火算法求解步骤：\par
Step1：初始化退火初始温度、终止温度、降温速度、初始解T\par
Step2：对T产生随机扰动，求解中心温度，判断是否满足约束条件，若不满足，重复该步骤\par
Step3：求解面积，若小于面积最小值，则更新面积最小值，并接受新解，否则，根据Metropolls准则接受新解\par
Step4：判断是否达到迭代次数L，若达到迭代次数，则降低温度，重置迭代次数\par
Step5：若温度小于终止温度，则输出结果，否则重复Step2、Step3、Step4、Step5

\begin{figure}[H]
	\centering
	\includegraphics[width=11.5cm]{问题三流程图.png} %[图片大小]{图片路径}
	\caption{算法流程图} %图片标题
\end{figure}


经过7447次迭代，求出最优参数
\begin{center}
	\begin{tabular}{|c|c|c|c|c|c|c|c|}
		\hline $T_0$& $T_1$ &$T_2$&$T_3$&$T_4$&$T_5$ & $v$ \\
		\hline 25&185&185&235&265&25&98\\
		\hline
	\end{tabular}
\end{center}


\begin{figure}[H]
	\centering
	\includegraphics[width=12cm]{untitled.png} %[图片大小]{图片路径}
	\caption{最小面积变化图} %图片标题
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=12cm]{第三问结果.png} %[图片大小]{图片路径}
	\caption{最优炉温曲线} %图片标题
\end{figure}
且其满足制程界限如下
\begin{table}[H]
	\centering
	\caption{制程界限}
	\begin{tabular}{ccc}
		\toprule[1.5pt]
		\makebox[0.3\textwidth][c]{参数}	&  \makebox[0.4\textwidth][c]{数值} \\ 
		\midrule[1pt]		
		名称                & 值      & 单位    \\ 
		\hline
		最大温度上升斜率         & 2.002  & \SI{}{\degreeCelsius}/s \\
		\hline
		最小温度上升斜率           & -1.770 & \SI{}{\degreeCelsius}/s \\
		\hline
		温度上升过程中在150\SI{}{\degreeCelsius}-190\SI{}{\degreeCelsius}的时间 & 75.5   &    s    \\
		\hline
		温度大于\SI{217}{\degreeCelsius}的时间	 &   67   & s  \\
		\hline
		峰值温度	   &   240.67   &      \SI{}{\degreeCelsius}\\
		\bottomrule[1.5pt]		
	\end{tabular}
\end{table}


% T0	T1	T2	T3	T4	T5	v
% 25	185	185	235	265	25	95

并给出相应的面积$Q_{\eta min}=525.44$。

\section{问题四的模型建立及求解}
\subsection{模型建立}
\subsubsection{定义指标}
对于以峰值温度为中心线的两侧超过217ºC的炉温曲线，关于中心线两侧，以时间间隔$\Delta t$划分。$t_{217}$表示上升过程中温度首次到达217摄氏度时的时间，$t_m$表示温度到达最大值时的时间，$t_i$表示温度介于217和最大值之间时的时间。我们定义两个指标 一个指标为$\psi$
\begin{equation}
	\begin{aligned}
		\begin{split}
		\left \{
			\begin{array}{ll}
		&\psi=\sum\left(u(i)-u\left(2 t_{m}-t_{i}\right)\right)^{2}\\
		&t_{i}=0.5 k+t_{217}\\
		&k \in N\\
		&t_{i} \leq t_{m}
			\end{array}
		\right.
		\end{split}
	\end{aligned}
\end{equation}
另一个指标为:
我们将第三问中的最优炉温曲线对应的阴影面积最小值$Q_{\eta min}$引入，定义一个指标$\zeta$
\begin{equation}
	\zeta=\frac{Q_\eta-Q_{\eta\mathrm{min}}}{Q_{\eta\mathrm{min}}}
\end{equation}
我们给定$\zeta\le0.2$的情况下，可以认为是最优炉温曲线


的情况下，可以认为是最优炉温曲线
\subsubsection{优化模型}
$\bullet$ 决策变量
\begin{equation}
	T=(T_1,T_2,T_3,T_4,v)
\end{equation}
 

$\bullet$ 目标函数
\begin{equation}
	Min \quad \Psi(T)
\end{equation}

$\bullet$ 约束条件
\begin{equation}
	\begin{aligned}
		\begin{split}
		\left \{
		\begin{array}{ll}
		    T_1\in[165,185]\\
			T_2\in[185,205]\\
			T_3\in[225,245]\\
			T_4\in[245,265]\\
			v\in[65,100]\\
			\frac{du}{dt}\in[-3,3]\\
			60 \leq\left. t\right|_{150 \leq u \leq 190} \leq 120 \\
			40 \leq\left. t\right|_{u \geq 270} \leq 90\\
			u_{max}\in[240,250]\\
			$\zeta\le 0.2$\\
		\end{array}
		\right.
		\end{split}
	\end{aligned}
\end{equation}

\subsection{问题求解}
使用模拟退火算法求解步骤
\par
Step1：初始化退火初始温度、终止温度、降温速度、初始解$T$
\par
Step2：对$T$产生随机扰动，求解中心温度，判断是否满足约束条件，若不满足，重复该步骤
\par
Step3：求解$\Psi$，若小于$\Psi $最小值，则更新$\Psi $最小值，并接受新解，否则，根据$Metropolls$准则接受新解
\par
Step4：判断是否达到迭代次数$L$，若达到迭代次数，则降低温度，重置迭代次数
\par
Step5：若温度小于终止温度，则输出结果，否则重复Step2、Step3、Step5
\par
求解出最优参数为(185,193,230,265,94)，面积527.57，\psi = 8010.31
\par
炉温曲线
\begin{figure}[H]
	\centering
	\includegraphics[width=14cm]{问题四温度曲线.png} 
	\caption{温度曲线} %图片标题
\end{figure}


\begin{center}
	\begin{array}{ccc}
		\hline \text { 名称 } & \text { 值 } & \text { 单位 } \\
		\hline \text { 最大温度上升斜率 } & 2.19 & { }^{\circ} \mathrm{C} / \mathrm{s} \\
		\hline \text { 最小温度下降斜率 } & -1.77 & { }^{\mathrm{o}} \mathrm{C} / \mathrm{s} \\
		\hline \text { 温度上升过程中在 } 150^{\circ} \mathrm{C} \sim 190^{\circ} \mathrm{C} \text { 的时间 } & 65.5 & \mathrm{~s} \\
		\hline \text { 温度大于 } 217^{\circ} \mathrm{C} \text { 的时间 } & 56.5 & \mathrm{~s} \\
		\hline \text { 峰值温度 } & 240 & { }^{\circ} \mathrm{C} \\
		\hline
	\end{array}
\end{center}


\section{模型检验}
\subsection{灵敏度分析}
在问题二当中，假设速度v在±5\%范围内波动，观察5个约束变量的变化，对其进行灵敏分析，得到如下所示的关系图：

\begin{figure}[H]
	\centering
	\includegraphics[width=14cm]{灵敏度分析.png} %[图片大小]{图片路径}
	\caption{灵敏度分析} %图片标题
\end{figure}
蓝色虚线线为每个约束条件的上下限，红色实线为v与约束变量的函数曲线，通过观察发现，当速度v在±5\%范围内波动时，温度变化的最大最小斜率均无明显变化，温度上升过程中在150ºC~190ºC的时间、温度大于217ºC的时间都随着速度的增大而减小，但都在上下限范围之内，峰值温度与速度v的函数关系图中，当v大于80.7时，峰值温度小于温度下限，不再满足约束条件
\subsection{模型优缺点}
\subsubsection{模型优点}
1.	在求解炉温曲线过程中，使用热传导方程、牛顿冷却定律，并分段拟合未知参数，能够较好的反应印刷电路板温度变化过程\par
2.	在求解问题三、问题四过程中，采用模拟退火方法，减少遍历次数，更快的寻找全局最优解
\subsubsection{模型缺点}
1.	未考虑回焊炉内各小温区温度对相邻小温度温度的影响
\par
2.	问题一中分段拟合求解未知参数，在粗略调整参数范围时，可能此范围含有局部最优解，而不含全局最优解
\begin{thebibliography}{1}

	\bibitem{}
	杨世铭,陶文铨.传热学(第四版)[M].北京:高等教育出版社, 2006.
	\bibitem{}
	李岩,赵立博,张伟,邱兆义.对回流焊炉温度设定的分析与优化[J].船电技术,2010,30(07):44-46.
\end{thebibliography}
% \bibliography{book}
% \bibliographystyle{IEEEtran}
% \bibliography{reference}

\newpage
%附录
\appendix
\section{求解炉内温度源代码}
\begin{lstlisting}[language=matlab]
function [tem] = stove_tem(x)
	%     T0 = 25;
	%     T1 = 175;
	%     T2 = 195;
	%     T3 = 235;
	%     T4 = 255;
	%     T5 = 25;
    
	T0 = 25.0;
	T1 = 173.0;
	T2 = 198.0;
	T3 = 230.0;
	T4 = 257.0;
	T5 = 25.0;

    if x >= 0 && x < 25
        tem = (T1 - T0) / 25 * (x - 0) + T0;
    elseif x >= 25 && x < 197.5
        tem = T1;
    elseif x >= 197.5 && x < 202.5
        tem = (T2 - T1) / 5 * (x - 197.5) + T1;
    elseif x >= 202.5 && x < 233
        tem = T2;
    elseif x >= 233 && x < 238
        tem = (T3 - T2) / 5 * (x - 233) + T2;
    elseif x >= 238 && x < 268.5 
        tem = T3;
    elseif x >= 268.5 && x < 273.5
        tem = (T4 - T3) / 5 * (x - 268.5) + T3;
    elseif x >= 273.5 && x < 339.5
        tem = T4;
    elseif x >= 339.5 && x < 344.5 
        tem = (T5 - T4) / 5 * (x - 339.5) + T4;
    elseif x >= 344.5 && x < 410.5
        tem = T5;
    elseif x >= 410.5 && x <= 435.5
        tem = (T0 - T5) / 5 * (x - 410.5) + T5;
    end
end
\end{lstlisting}

 \section{问题一源代码}
\begin{lstlisting}[language=matlab]
	clc; clear;

	stove_length = 435.5;
	speed1 = 78.0/60;

	L1_lim = 202.5;
	L2_lim = 233.0;
	L3_lim = 268.5;
	L4_lim = 339.5;
	L5_lim = stove_length;

	tim1 = L1_lim / speed1;
	tim2 = (L2_lim - L1_lim) / speed1;
	tim3 = (L3_lim - L2_lim) / speed1;
	tim4 = (L4_lim - L3_lim) / speed1;
	tim5 = (L5_lim - L4_lim) / speed1;

	T0 = 25.0;
	T1 = 173.0;
	T2 = 198.0;
	T3 = 230.0;
	T4 = 257.0;
	T5 = 25.0;

	Tim = stove_length / speed1;
	H = 0.15;
	dh = 0.01; 
	% dt = 0.00001;
	dt = 0.0001; 
	H_sum = floor(H / dh);
	T_sum = floor(Tim / dt);

	% lamda = 0.00003;
	lamda = 0.00002;

	u = zeros(H_sum, T_sum);

	lamda = 0.0000380;
	k_newton = 5000;
	lamda = lamda * (dt / (dh * dh));

	for i = 1 : 1 : H_sum
		u(i,1) = 25;
	end

	for k = 2 : 1 : T_sum
	%     u(1,k) = stove_tem(k * dt * speed1);
	%     u(H_sum,k) = stove_tem(k * dt * speed1);
		u(1,k) = -k_newton * dt * (u(1,k-1) - stove_tem(k * dt * speed1)) + u(1,k-1);
		u(H_sum,k) = -k_newton * dt * (u(H_sum,k-1) - stove_tem(k * dt * speed1)) + u(H_sum,k-1);
	end

	err = 0;
	k_cnt = 1;
	cnt = 1;
	for k = 2 : 1 : T_sum
		k_cnt = k_cnt + 1;
		for i = 2 : 1 : H_sum-1
			if k > tim1  / dt && k <= (tim1 +tim2) / dt
				lamda = 0.000045;
				lamda = lamda * (dt / (dh * dh));
			elseif k > (tim1 +tim2) / dt && k <= (tim1 +tim2 + tim3) / dt
				lamda = 0.000063;
				lamda = lamda * (dt / (dh * dh));
			elseif k > (tim1 +tim2 + tim3) / dt && k <= (tim1 +tim2+tim3+tim4)  / dt
				lamda = 0.000047;
				lamda = lamda * (dt / (dh * dh));
			elseif k > (tim1 +tim2+tim3+tim4)  / dt
				lamda = 0.000022;
				lamda = lamda * (dt / (dh * dh));
			end
			u(i,k) = lamda * u(i+1,k-1) + (1 - 2 * lamda) * u(i,k-1) + lamda * u(i-1,k-1);
			if u(i,k) >= 30 && i == floor(H_sum / 2) +1 && k_cnt == floor(0.5 / dt)
				pcb_calc_tem(cnt) = u(i,k);
				cnt = cnt+1;
			end
		end
		if k_cnt == floor(0.5 / dt)
			k_cnt = 0;
		end
	end         

	% plot(pcb_real_tem, 'r');

	plot(pcb_calc_tem, 'b');
	hold on;

	max_slop = 0;
	min_slop = 0;
	tim_up_150_190 = 0;
	tim_217 = 0;
	max_temp = 0;
	for i = 2 : 1 : length(pcb_calc_tem)
		slop = (pcb_calc_tem(i) - pcb_calc_tem(i - 1)) / 0.5;
		if slop > max_slop
			max_slop = slop;
		elseif slop < min_slop
			min_slop = slop
		end
		
		if pcb_calc_tem(i) > 217
			tim_217 = tim_217 + 0.5;
		end
		
		if pcb_calc_tem(i) > max_temp
			max_temp = pcb_calc_tem(i);
		end
		
		if pcb_calc_tem(i) >= 150 && pcb_calc_tem(i) <= 190 && pcb_calc_tem(i) > pcb_calc_tem(i-1)
			tim_up_150_190 = tim_up_150_190 + 0.5;
		end
	end

	% for i = 1 : 250
	%     line_lim(i,1) = tim1 * 2;
	%     line_lim(i,2) = (tim1 + tim2) * 2;
	%     line_lim(i,3) = (tim1 + tim2 + tim3) * 2;
	%     line_lim(i,4) = (tim1 + tim2 + tim3 + tim4) * 2;
	%     line_lim(i,5) = i;
	% end

	% plot(line_lim(:,1),line_lim(:,5),'--');
	% plot(line_lim(:,2),line_lim(:,5),'--');
	% plot(line_lim(:,3),line_lim(:,5),'--');
	% plot(line_lim(:,4),line_lim(:,5),'--');
 \end{lstlisting}
\section{问题二源代码}
\begin{lstlisting}[language=matlab]
	clc; clear;

	speed_cnt = 0;
	speed_length = 100 - 65 + 1;

	max_slop = zeros(speed_length,1);
	min_slop = zeros(speed_length,1);
	tim_up_150_190 = zeros(speed_length,1);
	tim_217 = zeros(speed_length,1);
	max_temp = zeros(speed_length,1);
	speed_ans = zeros(speed_length,1);
	for speed = 79 : 0.1 : 81
		speed_cnt = speed_cnt + 1;
		
		stove_length = 435.5;
		speed1 = speed / 60;

		L1_lim = 202.5;
		L2_lim = 233.0;
		L3_lim = 268.5;
		L4_lim = 339.5;
		L5_lim = stove_length;

		tim1 = L1_lim / speed1;
		tim2 = (L2_lim - L1_lim) / speed1;
		tim3 = (L3_lim - L2_lim) / speed1;
		tim4 = (L4_lim - L3_lim) / speed1;
		tim5 = (L5_lim - L4_lim) / speed1;

		T0 = 25.0;
		T1 = 182;
		T2 = 203;
		T3 = 237;
		T4 = 254;
		T5 = 25.0;

		Tim = stove_length / speed1;
		H = 0.15;
		dh = 0.01; 
		% dt = 0.00001;
		dt = 0.0001; 
		H_sum = floor(H / dh);
		T_sum = floor(Tim / dt);

		% lamda = 0.00003;
		lamda = 0.000038;

		U = zeros(H_sum, T_sum);

		k_newton = 5000;
		lamda = lamda * (dt / (dh * dh));

		for i = 1 : 1 : H_sum
			U(i,1) = 25;
		end
		
		for k = 2 : 1 : T_sum
			U(1,k) = -k_newton * dt * (U(1,k-1) - stove_tem(k * dt * speed1)) + U(1,k-1);
			U(H_sum,k) = -k_newton * dt * (U(H_sum,k-1) - stove_tem(k * dt * speed1)) + U(H_sum,k-1);
		end

		err = 0;
		k_cnt = 1;
		cnt = 1;
		for k = 2 : 1 : T_sum
			k_cnt = k_cnt + 1;
			for i = 2 : 1 : H_sum-1
				if k > tim1  / dt && k <= (tim1 +tim2) / dt
					lamda = 0.000045;
					lamda = lamda * (dt / (dh * dh));
				elseif k > (tim1 +tim2) / dt && k <= (tim1 +tim2 + tim3) / dt
					lamda = 0.000063;
					lamda = lamda * (dt / (dh * dh));
				elseif k > (tim1 +tim2 + tim3) / dt && k <= (tim1 +tim2+tim3+tim4)  / dt
					lamda = 0.000047;
					lamda = lamda * (dt / (dh * dh));
				elseif k > (tim1 +tim2+tim3+tim4)  / dt
					lamda = 0.000022;
					lamda = lamda * (dt / (dh * dh));
				end
				U(i,k) = lamda * U(i+1,k-1) + (1 - 2 * lamda) * U(i,k-1) + lamda * U(i-1,k-1);
				if U(i,k) >= 30 && i == floor(H_sum / 2) +1 && k_cnt == floor(0.5 / dt)
					pcb_calc_tem(cnt) = U(i,k);
					cnt = cnt+1;
				end
			end
			if k_cnt == floor(0.5 / dt)
				k_cnt = 0;
			end
		end         

		plot(pcb_calc_tem);
		hold on;


		for i = 2 : 1 : cnt - 1
			slop = (pcb_calc_tem(i) - pcb_calc_tem(i - 1)) / 0.5;
			if slop > max_slop(speed_cnt)
				max_slop(speed_cnt) = slop;
			elseif slop < min_slop(speed_cnt)
				min_slop(speed_cnt) = slop;
			end

			if pcb_calc_tem(i) > 217
				tim_217(speed_cnt) = tim_217(speed_cnt) + 0.5;
			end

			if pcb_calc_tem(i) > max_temp(speed_cnt)
				max_temp(speed_cnt) = pcb_calc_tem(i);
			end

			if pcb_calc_tem(i) >= 150 && pcb_calc_tem(i) <= 190 && pcb_calc_tem(i) > pcb_calc_tem(i-1)
				tim_up_150_190(speed_cnt) = tim_up_150_190(speed_cnt) + 0.5;
			end
		end
		
		if max_slop(speed_cnt) <= 3 && min_slop(speed_cnt) >= -3 ...
				&& tim_217(speed_cnt) >= 40 && tim_217(speed_cnt) <= 90 ...
				&& max_temp(speed_cnt) >= 240 && max_temp(speed_cnt) <= 250 ...
				&& tim_up_150_190(speed_cnt) >= 60 && tim_up_150_190(speed_cnt) <= 120
			speed_ans(speed_cnt) = speed;
		end
		
	end
	% for i = 1 : 250
	%     line_lim(i,1) = tim1 * 2;
	%     line_lim(i,2) = (tim1 + tim2) * 2;
	%     line_lim(i,3) = (tim1 + tim2 + tim3) * 2;
	%     line_lim(i,4) = (tim1 + tim2 + tim3 + tim4) * 2;
	%     line_lim(i,5) = i;
	% end

	% plot(line_lim(:,1),line_lim(:,5),'--');
	% plot(line_lim(:,2),line_lim(:,5),'--');
	% plot(line_lim(:,3),line_lim(:,5),'--');
	% plot(line_lim(:,4),line_lim(:,5),'--');
\end{lstlisting}

\section{问题三源代码}
\begin{lstlisting}[language=matlab]
	clc; clear;

	global T0;
	global T1;
	global T2;
	global T3;
	global T4;
	global T5;

	T0 = 25;
	T1 = 175;
	T2 = 195;
	T3 = 235;
	T4 = 255;
	T5 = 25;
	speed = 83;

	T0_last = T0;
	T1_last = T1;
	T2_last = T2;
	T3_last = T3;
	T4_last = T4;
	T5_last = T5;
	speed_last = speed;

	pcb_calc_tem = zeros(700,1);
	integral_all = zeros(1000,1);
	integral_all_min = zeros(1000,1);

	r = 0.95;
	L = 100;
	T_tui = 500;
	T_tui_min = 1;
	T_tui_cnt = 0;
	integral_now = 0;
	integral_last = 999999999999999;
	integral_max = 0;
	integral_min = 999999999999999999;
	dddd = 0;
	while T_tui > T_tui_min
		l_cnt = 1;
		while l_cnt < L
			dddd = dddd + 1
			T_tui_cnt
			% T_tui_cnt = T_tui_cnt + 1;
			% integral_last = integral_now;
			integral_now = 0;

			T1 = T1 + ceil(rand * 10) - 5;
			if T1 > 175 + 10
				T1 = 175 + 10;
			elseif T1 < 175 - 10
				T1 = 175 - 10;
			end

			T2 = T2 + ceil(rand * 10) - 5;
			if T2 > 195 + 10
				T2 = 195 + 10;
			elseif T2 < 195 - 10
				T2 = 195 - 10;
			end

			T3 = T3 + ceil(rand * 10) - 5;
			if T3 > 235 + 10
				T3 = 235 + 10;
			elseif T3 < 235 - 10
				T3 = 235 - 10;
			end

			T4 = T4 + ceil(rand * 10) - 5;
			if T4 > 255 + 10
				T4 = 255 + 10;
			elseif T4 < 255 - 10
				T4 = 255 - 10;
			end

			speed = speed + ceil(rand * 14) - 7;
			if speed > 100
				speed = 100;
			elseif speed < 65
				speed = 65;
			end


			stove_length = 435.5;
			speed1 = speed / 60;

			L1_lim = 202.5;
			L2_lim = 233.0;
			L3_lim = 268.5;
			L4_lim = 339.5;
			L5_lim = stove_length;

			tim1 = L1_lim / speed1;
			tim2 = (L2_lim - L1_lim) / speed1;
			tim3 = (L3_lim - L2_lim) / speed1;
			tim4 = (L4_lim - L3_lim) / speed1;
			tim5 = (L5_lim - L4_lim) / speed1;

			Tim = stove_length / speed1;
			H = 0.15;
			dh = 0.01; 
			% dt = 0.00001;
			dt = 0.0001; 
			H_sum = floor(H / dh);
			T_sum = floor(Tim / dt);

			u = zeros(H_sum, T_sum);

			lamda = 0.0000380;
			k_newton = 5000;
			lamda = lamda * (dt / (dh * dh));

			for i = 1 : 1 : H_sum
				u(i,1) = 25;
			end

			for k = 2 : 1 : T_sum
				x = k * dt * speed1;

				if x >= 0 && x < 25
					tem = (T1 - T0) / 25 * (x - 0) + T0;
				elseif x >= 25 && x < 197.5
					tem = T1;
				elseif x >= 197.5 && x < 202.5
					tem = (T2 - T1) / 5 * (x - 197.5) + T1;
				elseif x >= 202.5 && x < 233
					tem = T2;
				elseif x >= 233 && x < 238
					tem = (T3 - T2) / 5 * (x - 233) + T2;
				elseif x >= 238 && x < 268.5 
					tem = T3;
				elseif x >= 268.5 && x < 273.5
					tem = (T4 - T3) / 5 * (x - 268.5) + T3;
				elseif x >= 273.5 && x < 339.5
					tem = T4;
				elseif x >= 339.5 && x < 344.5 
					tem = (T5 - T4) / 5 * (x - 339.5) + T4;
				elseif x >= 344.5 && x < 410.5
					tem = T5;
				elseif x >= 410.5 && x <= 435.5
					tem = (T0 - T5) / 5 * (x - 410.5) + T5;
				else
					tem = 0;
				end

				u(1,k) = -k_newton * dt * (u(1,k-1) - tem) + u(1,k-1);
				u(H_sum,k) = -k_newton * dt * (u(H_sum,k-1) - tem) + u(H_sum,k-1);
			end

			err = 0;
			k_cnt = 1;
			cnt = 1;
			for k = 2 : 1 : T_sum
				k_cnt = k_cnt + 1;
				for i = 2 : 1 : H_sum-1
					if k > tim1  / dt && k <= (tim1 +tim2) / dt
						lamda = 0.000045;
						lamda = lamda * (dt / (dh * dh));
					elseif k > (tim1 +tim2) / dt && k <= (tim1 +tim2 + tim3) / dt
						lamda = 0.000063;
						lamda = lamda * (dt / (dh * dh));
					elseif k > (tim1 +tim2 + tim3) / dt && k <= (tim1 +tim2+tim3+tim4)  / dt
						lamda = 0.000047;
						lamda = lamda * (dt / (dh * dh));
					elseif k > (tim1 +tim2+tim3+tim4)  / dt
						lamda = 0.000022;
						lamda = lamda * (dt / (dh * dh));
					end
					u(i,k) = lamda * u(i+1,k-1) + (1 - 2 * lamda) * u(i,k-1) + lamda * u(i-1,k-1);
					if u(i,k) >= 30 && i == floor(H_sum / 2) +1 && k_cnt == floor(0.5 / dt)
						pcb_calc_tem(cnt) = u(i,k);
						cnt = cnt+1;
					end
				end
				if k_cnt == floor(0.5 / dt)
					k_cnt = 0;
				end
			end         
		%     plot(pcb_calc_tem, 'b');
		%     hold on;

			max_slop = 0;
			min_slop = 0;
			tim_up_150_190 = 0;
			tim_217 = 0;
			max_temp = 0;
			for i = 2 : 1 : cnt - 1
				slop = (pcb_calc_tem(i) - pcb_calc_tem(i - 1)) / 0.5;
				if slop > max_slop
					max_slop = slop;
				elseif slop < min_slop
					min_slop = slop;
				end

				if pcb_calc_tem(i) > 217
					tim_217 = tim_217 + 0.5;
				end

				if pcb_calc_tem(i) > max_temp
					max_temp = pcb_calc_tem(i);
				end

				if pcb_calc_tem(i) >= 150 && pcb_calc_tem(i) <= 190 && pcb_calc_tem(i) > pcb_calc_tem(i-1)
					tim_up_150_190 = tim_up_150_190 + 0.5;
				end

				if(pcb_calc_tem(i) >= 217)
					integral_now = integral_now + 2 / pi * 0.5 * (pcb_calc_tem(i) - 217);
				end
			end

			if max_slop <= 3 && min_slop >= -3 ...
					&& tim_217 >= 40 && tim_217 <= 90 ...
					&& max_temp >= 240 && max_temp <= 250 ...
					&& tim_up_150_190 >= 60 && tim_up_150_190 <= 120
					l_cnt = l_cnt + 1;
			else
				continue;
			end    

			T_tui_cnt = T_tui_cnt + 1
			integral_all(T_tui_cnt) = integral_now;

			dE = integral_now - integral_last
			if dE <= 0
				if integral_min > integral_now
					integral_min = integral_now;
					t0_ans = T0;
					t1_ans = T1;
					t2_ans = T2;
					t3_ans = T3;
					t4_ans = T4;
					t5_ans = T5;
					speed_ans = speed;
				end
			else
				if exp(-dE / T_tui) > rand

				else
					T0 = T0_last;
					T1 = T1_last;
					T2 = T2_last;
					T3 = T3_last;
					T4 = T4_last;
					T5 = T5_last;
					speed = speed_last;
				end
			end

			integral_all_min(T_tui_cnt) = integral_min;

			integral_last = integral_now;

			T0_last = T0;
			T1_last = T1;
			T2_last = T2;
			T3_last = T3;
			T4_last = T4;
			T5_last = T5;
			speed_last = speed;
		end
		
		T_tui = r * T_tui;
		% for i = 1 : 250
		%     line_lim(i,1) = tim1 * 2;
		%     line_lim(i,2) = (tim1 + tim2) * 2;
		%     line_lim(i,3) = (tim1 + tim2 + tim3) * 2;
		%     line_lim(i,4) = (tim1 + tim2 + tim3 + tim4) * 2;
		%     line_lim(i,5) = i;
		% end

		% plot(line_lim(:,1),line_lim(:,5),'--');
		% plot(line_lim(:,2),line_lim(:,5),'--');
		% plot(line_lim(:,3),line_lim(:,5),'--');
		% plot(line_lim(:,4),line_lim(:,5),'--');
	end
\end{lstlisting}

\section{问题四源代码}
\begin{lstlisting}[language=matlab]
	clc; clear;

	global T0;
	global T1;
	global T2;
	global T3;
	global T4;
	global T5;

	T0 = 25;
	T1 = 175;
	T2 = 195;
	T3 = 235;
	T4 = 255;
	T5 = 25;
	speed = 83;

	T0_last = T0;
	T1_last = T1;
	T2_last = T2;
	T3_last = T3;
	T4_last = T4;
	T5_last = T5;
	speed_last = speed;

	r = 0.95;
	L = 10;
	T_tui = 1000;
	T_tui_min = 1;
	T_tui_cnt = 0;
	integral_now = 0;
	integral_last = 999999999999999;
	integral_max = 0;
	integral_min = 999999999999999999;
	err_temp_last = 99999999999999999;
	err_temp_min = 99999999999999999;
	dddd = 0;
	while T_tui > T_tui_min
		l_cnt = 1;
		while l_cnt < L
			dddd = dddd + 1
			T_tui_cnt
			% T_tui_cnt = T_tui_cnt + 1;
			% integral_last = integral_now;
			integral_now = 0;

			T1 = T1 + ceil(rand * 10) - 5;
			if T1 > 175 + 10
				T1 = 175 + 10;
			elseif T1 < 175 - 10
				T1 = 175 - 10;
			end

			T2 = T2 + ceil(rand * 10) - 5;
			if T2 > 195 + 10
				T2 = 195 + 10;
			elseif T2 < 195 - 10
				T2 = 195 - 10;
			end

			T3 = T3 + ceil(rand * 10) - 5;
			if T3 > 235 + 10
				T3 = 235 + 10;
			elseif T3 < 235 - 10
				T3 = 235 - 10;
			end

			T4 = T4 + ceil(rand * 10) - 5;
			if T4 > 255 + 10
				T4 = 255 + 10;
			elseif T4 < 255 - 10
				T4 = 255 - 10;
			end

			speed = speed + ceil(rand * 14) - 7;
			if speed > 100
				speed = 100;
			elseif speed < 65
				speed = 65;
			end


			stove_length = 435.5;
			speed1 = speed / 60;

			L1_lim = 202.5;
			L2_lim = 233.0;
			L3_lim = 268.5;
			L4_lim = 339.5;
			L5_lim = stove_length;

			tim1 = L1_lim / speed1;
			tim2 = (L2_lim - L1_lim) / speed1;
			tim3 = (L3_lim - L2_lim) / speed1;
			tim4 = (L4_lim - L3_lim) / speed1;
			tim5 = (L5_lim - L4_lim) / speed1;

			Tim = stove_length / speed1;
			H = 0.15;
			dh = 0.01; 
			% dt = 0.00001;
			dt = 0.0001; 
			H_sum = floor(H / dh);
			T_sum = floor(Tim / dt);

			u = zeros(H_sum, T_sum);

			lamda = 0.0000380;
			k_newton = 5000;
			lamda = lamda * (dt / (dh * dh));

			for i = 1 : 1 : H_sum
				u(i,1) = 25;
			end

			for k = 2 : 1 : T_sum
				x = k * dt * speed1;

				if x >= 0 && x < 25
					tem = (T1 - T0) / 25 * (x - 0) + T0;
				elseif x >= 25 && x < 197.5
					tem = T1;
				elseif x >= 197.5 && x < 202.5
					tem = (T2 - T1) / 5 * (x - 197.5) + T1;
				elseif x >= 202.5 && x < 233
					tem = T2;
				elseif x >= 233 && x < 238
					tem = (T3 - T2) / 5 * (x - 233) + T2;
				elseif x >= 238 && x < 268.5 
					tem = T3;
				elseif x >= 268.5 && x < 273.5
					tem = (T4 - T3) / 5 * (x - 268.5) + T3;
				elseif x >= 273.5 && x < 339.5
					tem = T4;
				elseif x >= 339.5 && x < 344.5 
					tem = (T5 - T4) / 5 * (x - 339.5) + T4;
				elseif x >= 344.5 && x < 410.5
					tem = T5;
				elseif x >= 410.5 && x <= 435.5
					tem = (T0 - T5) / 5 * (x - 410.5) + T5;
				else
					tem = 0;
				end

				u(1,k) = -k_newton * dt * (u(1,k-1) - tem) + u(1,k-1);
				u(H_sum,k) = -k_newton * dt * (u(H_sum,k-1) - tem) + u(H_sum,k-1);
			end

			err = 0;
			k_cnt = 1;
			cnt = 1;
			for k = 2 : 1 : T_sum
				k_cnt = k_cnt + 1;
				for i = 2 : 1 : H_sum-1
					if k > tim1  / dt && k <= (tim1 +tim2) / dt
						lamda = 0.000045;
						lamda = lamda * (dt / (dh * dh));
					elseif k > (tim1 +tim2) / dt && k <= (tim1 +tim2 + tim3) / dt
						lamda = 0.000063;
						lamda = lamda * (dt / (dh * dh));
					elseif k > (tim1 +tim2 + tim3) / dt && k <= (tim1 +tim2+tim3+tim4)  / dt
						lamda = 0.000047;
						lamda = lamda * (dt / (dh * dh));
					elseif k > (tim1 +tim2+tim3+tim4)  / dt
						lamda = 0.000022;
						lamda = lamda * (dt / (dh * dh));
					end
					u(i,k) = lamda * u(i+1,k-1) + (1 - 2 * lamda) * u(i,k-1) + lamda * u(i-1,k-1);
					if u(i,k) >= 30 && i == floor(H_sum / 2) +1 && k_cnt == floor(0.5 / dt)
						pcb_calc_tem(cnt) = u(i,k);
						cnt = cnt+1;
					end
				end
				if k_cnt == floor(0.5 / dt)
					k_cnt = 0;
				end
			end         
		%     plot(pcb_calc_tem, 'b');
		%     hold on;

			max_slop = 0;
			min_slop = 0;
			tim_up_150_190 = 0;
			tim_217 = 0;
			tim_up_217_index = 0;
			max_temp = 0;
			max_temp_index = 0;
			for i = 2 : 1 : cnt - 1
				slop = (pcb_calc_tem(i) - pcb_calc_tem(i - 1)) / 0.5;
				if slop > max_slop
					max_slop = slop;
				elseif slop < min_slop
					min_slop = slop;
				end

				if pcb_calc_tem(i) > 217
					tim_217 = tim_217 + 0.5;
				end

				if pcb_calc_tem(i) >= 217 && pcb_calc_tem(i-1) < 217
					tim_up_217_index = i;
				end

				if pcb_calc_tem(i) > max_temp
					max_temp = pcb_calc_tem(i);
					max_temp_index = i;
				end

				if pcb_calc_tem(i) >= 150 && pcb_calc_tem(i) <= 190 && pcb_calc_tem(i) > pcb_calc_tem(i-1)
					tim_up_150_190 = tim_up_150_190 + 0.5;
				end

				if(pcb_calc_tem(i) >= 217)
					integral_now = integral_now + 2 / pi * 0.5 * (pcb_calc_tem(i) - 217);
				end
			end
			integral_now
			if max_slop <= 3 && min_slop >= -3 ...
					&& tim_217 >= 40 && tim_217 <= 90 ...
					&& max_temp >= 240 && max_temp <= 250 ...
					&& tim_up_150_190 >= 60 && tim_up_150_190 <= 120 ...
					&& (integral_now - 525.44) / 525.44 < 0.2 && (integral_now - 525.44) / 525.44 > -0.2
				l_cnt = l_cnt + 1;
			else
				continue;
			end    

			err_temp = 0;
			for i = tim_up_217_index : 1 : max_temp_index - 1
				err_temp = err_temp + (pcb_calc_tem(i) - pcb_calc_tem(2 * max_temp_index - i))^2;
			end

			T_tui_cnt = T_tui_cnt + 1
			
			err_temp_all(T_tui_cnt) = err_temp;
			dE = err_temp - err_temp_last
			if dE <= 0
				if err_temp_min > err_temp
					err_temp_min = err_temp;
					t0_ans = T0;
					t1_ans = T1;
					t2_ans = T2;
					t3_ans = T3;
					t4_ans = T4;
					t5_ans = T5;
					speed_ans = speed;
				end
			else
				if exp(-dE / T_tui) > rand

				else
					T0 = T0_last;
					T1 = T1_last;
					T2 = T2_last;
					T3 = T3_last;
					T4 = T4_last;
					T5 = T5_last;
					speed = speed_last;
				end
			end
			err_temp_min_all(T_tui_cnt) = err_temp_min;
			integral_last = integral_now;

			T0_last = T0;
			T1_last = T1;
			T2_last = T2;
			T3_last = T3;
			T4_last = T4;
			T5_last = T5;
			speed_last = speed;
		end
		T_tui = r * T_tui;
		% for i = 1 : 250
		%     line_lim(i,1) = tim1 * 2;
		%     line_lim(i,2) = (tim1 + tim2) * 2;
		%     line_lim(i,3) = (tim1 + tim2 + tim3) * 2;
		%     line_lim(i,4) = (tim1 + tim2 + tim3 + tim4) * 2;
		%     line_lim(i,5) = i;
		% end

		% plot(line_lim(:,1),line_lim(:,5),'--');
		% plot(line_lim(:,2),line_lim(:,5),'--');
		% plot(line_lim(:,3),line_lim(:,5),'--');
		% plot(line_lim(:,4),line_lim(:,5),'--');
	end
\end{lstlisting}
